Finite Difference Time Domain (FDTD) Simulations of Electromagnetic Wave 

Propagation Using a Spreadsheet 
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We describe a simple and intuitive implementation of the method of finite difference time domain 
simulations for propagating electromagnetic waves using the simplest possible tools available in Mi- 
crosoft Excel. The method overcomes the usual obstacles of familiarity with programming languages 
as it relies on little more than the cut and paste features that are standard in Excel. Avenues of 
exploration by students are proposed and sample graphs are included. The pedagogical effectiveness 
of the implementation was tested during an Independent Activities Period class, composed of 80'/, 
freshmen, at MIT, and yielded positive results. 



I. INTRODUCTION 

Here we outline a simple demonstrative example of the 
method of finite difference time domain (FDTD) sim- 
ulations of electromagnetic wave propagation appropri- 
ate for undergraduates in an introductory electricity and 
magnetism course or advanced high school science stu- 
dents using Microsoft Excel, a robust software package 
that offers advanced graphics, numeric, and animation 
capabilities requiring minimal computer experienceiS*^ 
Only two rows in the spreadsheet, one for initial elec- 
tric field and another for initial magnetic field, need be 
specified in order to initiate the simulation, and graphs 
of the spatially and temporally evolving waveforms can 
be updated in real time as the simulation is carried out 
in Excel. We begin with a review of the basic tenets of 
the FDTD method, proceed into an outline for the im- 
plementation in Microsoft Excel, and conclude with some 
suggested exercises. 



II. METHODS 

A. FDTD basics 

The starting point for an FDTD simulation are 
Maxwell's equations, which are repeated here for the case 
of one dimensional free space propagation in time {t) and 
space (z) with no sources or sinks for magnetic or electric 
fields B or E respectively for the corresponding material 
responses H or D, 



dt 



dHy 
dt 



1 dHy 

Co dz 

1 dE^ 
liQ dz 



(1) 



(2) 



V ■E = 0, 



V • B = 0. 



If we choose the x-direction for the polarization of the 
electric field and the z-direction for direction of propaga- 
tion, then it follows that the magnetic field is y-polarized 
as indicated in the subscripts on the fields above. To dis- 
cretize the equations of propagation, central difference 
approximations are taken for the derivatives of time and 
space. Temporal steps are indexed by the integer n and 
related to continuous time by the relation t = nAt, and 
spatial steps are indexed by the integer k and related to 
continuous space by the relation z = kAz. The temporal 
discretization method using central differences can then 
be designated for variable X as 



dX 
'dt 



AX 
'At 



X{n + l/2) - X{n-l/2) 



At 



and spatial discretization as 



dX 

dz 



AX 
~A^ 



X{k + l/2) - X{k-l/2) 



(5) 



(6) 



In FDTD, we need only consider the two curl equations 
Hand 121 above because the divergence conditions 01 and ^ 
can be satisfied implicitly by interleaving the electric and 
magnetic field components in space in what has come to 
be known as a Yee cell.^ A consequence of the spatial 
interleaving is that the fields must also be interleaved in 
time, known as "leapfrog", since the temporal response 
of one field is proportional to the spatial variation of the 
other at the previous time step. Implementing eqs. ^ 
and as indicated in eqs. 10 and © and interleaving 
spatially and temporally yields the difference equation 
version of Maxwell's equations: 
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At a given position in space, the field at each current 
time step can be calculated from the previous values of 
the fields. Solving for the latest field at the latest time 
step in eqs. and ^ yields: 



Err 



At 
eoAz 



Hy{k 



(9) 
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In this fashion, known as Euler forward, the solution 
of Maxwell's equations proceeds much the same way that 
we envision electromagnetic wave propagation-an electric 
field induces a magnetic field, which induces and electric 
field, ad infinitum. The time step and grid size parame- 
ters are chosen based on the propagating wave frequency 
and wavelength. For stability, the general rule is that at 
least ten grid points sum to less than the smallest wave- 
length \min considered, and the Courant condition then 
determines the time stept^ 



Az < 



10 



At < 



Az 
2^' 



(11) 



(12) 



where cq is the speed of light in vacuum. 

The computer algorithm implementation of this pro- 
cedure, eqs. (jni and (fTn|l . requires two one-dimensional 
arrays in the spatial coordinate, one for Ex and one for 
Hy, to be allocated in memory. After imposing an initial 
condition, in the form of an initial electric or magnetic 
field, each time step is calculated by the following pre- 
scription: 



At 



t ~ tn + n- 



Ex[k] = Ex[k] + -{Hy[k - 1] - Hy[k]} 

At 

t = to + n — 



B. FDTD spreadsheet algorithm 

The FDTD methodology for one dimension can be im- 
plemented in spreadsheet format, here in Microsoft Ex- 
cel, using simple cell formulas and cut-and-paste fea- 
tures. This has an advantage over other pedagogical 
approaches^ in that no programming experience is re- 
quired. The starting point is the algorithm in where 
columns represent spatial steps and pairs of rows (one 
for E and one for H) represent time steps. The algo- 
rithm is illustrated graphically in the center section of 
figure 1. For any time step (any pair of rows of E and 
H), the mapping of is: Ex{column) / Hy{column). 
Since the electric field value from one prior cell is needed 
to compute the curl. Ex begins at column one and ends 
at column k + 1. Similarly, the magnetic field begins 
at column zero but ends at column fc, since its value at 
one subsequent cell is needed when computing the curl. 
The first two rows [t = 0) are for the initial conditions. 
The algorithm computes the remaining cells. The stu- 
dent need only type in the formulas for E and H for the 
second time step and the first two columns; cut and paste 
may be used to complete the remaining spatial columns 
and temporal rows. 

The procedure for implementing the algorithm is as 
follows: 

1. As a visual aid, enter the spatial coordinate (1, 2, 3, 
etc) in the topmost row of the spreadsheet starting 
at column B. After typing in the first two spatial 
coordinates, select cells B and C and drag the "fill" 
handle to encompass as large a problem space as 
needed. 

2. Type 'ExO' in cell A2 and 'HyO' in cell A3, where 
the number following the field component refers to 
the time step. Drag the "fill" handle down the col- 
umn, as in 1, to include the desired number of time 
steps. 

3. Highlight the leftmost and rightmost columns of 
the problem space. 

4. Fill in all of row 2 and 3 (time step 0) to the end 
of the problem space with zeros and highlight to 
indicate that this is the initial condition for the 
fields. 



Hy[k] = Hy[k] + -{Ex[k] ~Ex[k+ 1]} 



(13) 



where we have adopted the normalized fields [E — 
y/eoiMjE) to simplify the code, and included the stabil- 
ity conditions, eqs. (|ll|l and l|12|l . within the constant 
preceding the curli^ The algorithm (eqn. H13|l) is iter- 
ated for the desired number of time steps. Note that 
iteration of H13|l results in an implicit time formulation, 
i.e. time is not made explicit in the equations for the field 
and only appears in the algorithm above for bookkeeping 
purposes. 



5. In the leftmost column and first time step for Ex, 
type '=B2'. Do the same for the rightmost col- 
umn. This imposes a perfect electric conductor ra- 
diation boundary condition on the problem space, 
which means that impinging waves will reflect from 
the problem space boundaries as they would from 
a mirror. 

6. In the second spatial position of Ex for time step 
one (ceU C4), type the formula '=C2-K0.5*(B3-C3)' 
and press enter. 
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7. Select the cell from the previous step and drag the 
"fill" handle to the last un-highlighted spatial col- 
umn in the row. 

8. In the first spatial position of Hy for time step one 
(cell C5) type the formula '=B2+0.5*(B4-C4)' and 
press enter. These last two steps define the com- 
puting steps of the algorithm. 

9. Select the cell from the previous step, and drag 
the "fill" handle to the last un-highlighted spatial 
column in the row. 

10. Finally, select all the columns for time step 1 [E^ 
and Hy), and drag the "fill" handle through the last 
time step in the simulation. Test that it works by 
entering a 1 in the initial condition regions and see 
that all cells automatically update; note that the 
only valid regions for initialization are at the zero 
time step {E^ and/or Hy) and between the first and 
last columns of the problem space (highlighted). 

The FDTD Excel code will now update all fields in 
response to the initial conditions entered by the user. To 
graphically visualize the simulation output, a single row 
can be selected and graphed using the insert > chart > 
(xy) scatter menu item. This outputs the spatial field 
pattern at the moment in time designated by the row 
number. To graph temporal evolution at a single point in 
space, it is recommended to make a new field mesh with 
only one field component (either Ex or Hy) on a separate 
sheet and graph from that; otherwise, both the E and H 
fields will appear in the graphs. This can be done using 
cut-and-paste. An example spatial and temporal graph 
of two counter-propagating pulses is illustrated in figure 
3. 



III. APPLICATIONS OF THE ALGORITHM 

A. Power flow 

The direction of propagation of light is dictated by 
the phase between the electric and magnetic fields. For 
transverse electromagnetic (TEM) wave propagation the 
familiar right-hand rule applied by curling E into H indi- 
cates the direction of power fiow. Using our spreadsheet 
code, this can be demonstrated by initializing both the 
electric and magnetic field as a sinusoidal waveform mod- 
ulated by a Gaussian envelope. Introducing a phase shift 
of either or tt in the magnetic field allows the direction 
of propagation to be controlled. If the electric {E^) and 
magnetic {Hy) field are both positive at the same points 
in space and time, then the pulse travels in the positive 
z direction; if the signs are opposite then propagation is 
in the opposite direction as illustrated in figure 2. 



B. Guided wave optics 

Some of the most interesting applications of electrody- 
namics occur when constraints are placed on the fields. 
The perfect electric conductor radiation boundary condi- 
tions form the walls of a resonator when the wavelength 
in the problem space (resonator width) is of similar size 
to the problem space itself. With PEC walls, no light 
wave with a wavelength longer than twice the resonator 
width will propagate. This is referred to as a cut-off 
wavelength. Also, since the electric field must be zero at 
the boundary and only nodes in a plane wave may sat- 
isfy this, only wavelengths that are half-integer multiples 
of the resonator width can propagate. Resonator modes 
can be investigated by setting the initial conditions to a 
plane wave with nodes at the boundaries is illustrated in 
figure 3. 



C. Transmission and reflection 

By introducing a relative permittivity into equa- 
tion lO , which changes the index of refraction and hence 
the wave propagation speed, the spreadsheet code can 
illustrate reflection and transmission at a material inter- 
face. In terms of the normalized fields, the constant 1/2 
becomes l/2ey.- Small errors in the transmitted and re- 
flected fields are to be expected, but these will decrease as 
the number of mesh points per wavelength increases. Fig- 
ure 4 illustrates this process for an air (e^ = 1) /material 
{Sr = 2) interface. 



IV. FURTHER EXERCISES AND SUGGESTED 
ACTIVITIES 

The following exercises are recommended. The spread- 
sheet code for the examples is available online, but it is 
recommended that students augment the code on their 
own, because writing and tinkering with the code is the 
best way to understand how it works. Exercises are listed 
in order of increasing difficulty, or in a manner such that 
each exercise depends only on changes already adminis- 
tered. 



A. Fundamental 

1. Use trigonometric and exponential mathematical 
functions to generate the initial conditions and re- 
produce the results in figure 2. 

2. Study the cavity modes of a one-dimensional res- 
onator, figure 3, by introducing standing waves into 
the initial conditions. What happens when the 
wavelength becomes longer than the problem space 
size (the width of the resonator)? 
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3. Experiment with radiation boundary conditions. 
Change the radiation boundary condition from per- 
fect electric conductor to perfect magnetic conduc- 
tor (Hj^ is zero at the problem space boundaries). 
What is different about reflections from these two 
kinds of surfaces? 

4. Implement a periodic boundary condition for un- 
bounded propagation. 

5. Reproduce the reflection and transmission simula- 
tion in figure 4, and then change the interface from 
air /material to material/ air. What changes? 

B. Advanced 

1. Experiment with spatially periodic variation in 
the relative permittivity. This should produce 
frequency-dependent filtering, and for some condi- 
tions should produce " photonic bandgap" materials 
with no propagation in certain frequency regions. 

2. Introduce loss into the simulation by including non- 
zero conductivity in Maxwell's equations and im- 
plementing it into the FDTD code. 

3. Implement an absorbing boundary condition. 



Independent Activities Period in a short course^ that met 
for two hours on each of two consecutive days. This brief 
introduction was sufficient for students to gain substan- 
tial proficiency in elementary simulations like those il- 
lustrated above. The class enjoyed unusually high at- 
tendance, especially among students preparing for or re- 
cently involved in electromagnetism coursework. Instruc- 
tor (D.W.W.) observation and student feedback after the 
course indicated considerable satisfaction with insights 
and proficiency gained. Several student responses pin- 
pointed the graphical form of the spreadsheet based al- 
gorithm as the key merit of the approach. The benefit of 
keeping track of a small number of unit cells in order to 
see how the algorithm works is useful for beginning and 
intermediate students of computer simulation, numerical 
recipes, and electricity and magnetism. The automatic 
updating of Excel graphs makes this an effective learning 
tool for students interested in the basics of electromag- 
netic wave propagation, as the consequences of changes 
in initial conditions are illustrated graphically and in- 
stantly. Finally, the methodology learned in this rapid 
introduction may stimulate student interest in more ad- 
vanced FDTD simulation techniques and their broad re- 
search applications JiiS 



V. CONCLUSION 

We have formulated a simple implementation of the 
FDTD method of propagating electromagnetic waves us- 
ing basic features available in Microsoft Excel, and we 
have presented some illustrative examples. The imple- 
mentation was tested on freshmen at MIT during its 2004 
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FIG. hGraphjcal illustration of the Spneadsheet iTinplementatlon of F DTD. Center: Com ponents of 
the spreadsheet "eld^ an-d a g^aphic^l rriap of the "eld up<Jate dependencies. Exterior: Illustration 
of the I mplemen tation 5 p? 1 n th? le:< t 



FIG. 2: INustralion of the nelatronship betweer the phas.es of E and H and the direction of 
pfOp^g^rlori.a) Left-tr^vdlrty piAi^ gf^nersiecJ by Introducing d p\/2 phui^nUtfX m th^-irltJ^IH 
"eld wilh respect to E, b) Righl-lraveNng pulse generated by initializing E and H in pha^e With 
e^ch other. Both graphs ihow the "elds afterthe s^i™ (arbitral^) -nurnher of tinn* ste-ps after 
g^rteratlon.Thf^ ■irnall^r pgls^s In a) ^neH b) are rernrhants of irnperfe<l Irililal condirl^>ris.Th$ 
armw indicates the direction of propagation and the dotted line down the center demarks the. 
5tart In g i?o I nt of the pu Ise,. 



Lowest energy symmetric and anttsymmetric (bold) cavity mode for a 50 unit wide resonant 
cwky with perfect el^tnc conductor walls. Flgure^s d I splays the jpatial dtstributlon of ihe etectdc 
and magnecic "elds in the cavity after 200 lime steps. 
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